package ru.vsu.amm.algebra;

import java.util.Arrays;

/**
 * Матрица, её характеристики и связанные алгоритмы
 * @author Иванов Илья
 * @since 02 января 2013
 */
public class Matrix {
	
	private double[][] elems;
	
	private int rows;
	private int columns;
	
	private static final int EIGENVECTOR_MAX_ITERATION = 10000;
	
	public Matrix(int dim) {
		this(dim, dim);
	}
	
	public Matrix(int rows, int columns) {
		if (rows < 0 || columns < 0 || rows * columns == 0 && (rows > 0 || columns > 0))
			throw new IncorrectDimensionException("Некорректное значение размерности: " + rows + "x" + columns);
		
		init(this, rows, columns);
	}
	
	private Matrix() {}
	
	private static Matrix init(Matrix matrix, int rows, int columns) {
		matrix.elems = new double[rows][columns];
		matrix.rows = rows;
		matrix.columns = columns;
		
		return matrix;
	}
	
	public Matrix copy() {
		Matrix matrix = init(new Matrix(), rows(), columns());
		for (int i = 0; i < rows(); i++)
			for (int j = 0; j < columns(); j++)
				matrix.set(get(i, j), i, j);
		
		return matrix;
	}
	
	public void copyFrom(Matrix matrix) {
		if (matrix.rows() != rows() || matrix.columns() != columns())
			throw new IncompatibleDimensionsException();
		
		for (int i = 0; i < rows(); i++)
			for (int j = 0; j < columns(); j++)
				set(matrix.get(i, j), i, j);
	}
	
	public int rows() {
		return rows;
	}
	
	public int columns() {
		return columns;
	}
	
	public boolean isSquare() {
		return rows() == columns();
	}
	
	public void set(double value, int row, int column) {
		elems[row][column] = value;
	}
	
	public double get(int row, int column) {
		return elems[row][column];
	}
	
	/**
	 * @param row Индекс строки.
	 * @return Копия заданной строки.
	 */
	public double[] getRow(int row) {
		if (row < 0 || row >= rows())
			throw new IndexOutOfBoundsException("Индекс строки выходит за границы размерности матрицы");
		
		return Arrays.copyOf(elems[row], columns());
	}
	
	/**
	 * @param column Индекс столбца.
	 * @return Копия заданного столбца.
	 */
	public double[] getColumn(int column) {
		if (column < 0 || column >= columns())
			throw new IndexOutOfBoundsException("Индекс столбца выходит за границы размерности матрицы");
		
		double[] columnElems = new double[rows()];
		for (int i = 0; i < columnElems.length; i++)
			columnElems[i] = get(i, column);
		
		return columnElems;
	}
	
	public void zero() {
		for (int i = 0; i < rows(); i++)
			for (int j = 0; j < columns(); j++)
				set(0, i, j);
	}
	
	public void identity() {
		if (!isSquare())
			throw new IllegalOperationException("Единичная матрица должна быть квадратной");
		
		zero();
		for (int k = 0; k < rows(); k++)
			set(1, k, k);
	}
	
	public void add(Matrix matrix) {
		performElemOperation(this, ElemOperation.ADDITION, matrix, this, false);
	}
	
	public void subtract(Matrix matrix) {
		performElemOperation(this, ElemOperation.SUBTRACTION, matrix, this, false);
	}
		
	public void multiply(double alpha) {
		for (int i = 0; i < rows(); i++)
			for (int j = 0; j < columns(); j++)
				set(alpha * get(i, j), i, j);
	}
	
	public Matrix normalize() {
		if (rows() > 1 && columns() > 1)
			throw new IllegalOperationException("Нормировка определена только для векторов");
		
		if (rows() == 1) {
			double sum = 0;
			for (int j = 0; j < columns(); j++)
				sum += get(0, j);
			for (int j = 0; j < columns(); j++)
				set(get(0, j) / sum, 0, j);
		} else {
			double sum = 0;
			for (int i = 0; i < rows(); i++)
				sum += get(i, 0);
			for (int i = 0; i < rows(); i++)
				set(get(i, 0) / sum, i, 0);
		}
		
		return this;
	}
	
	/**
	 * @param precision Точность вычислений.
	 * @return Собственный вектор, отвечающий максимальному по модулю собственному значению.
	 */
	public Matrix maxEigenvector(double precision) {
		if (!isSquare())
			throw new IllegalOperationException("Поиск собственных векторов определён только для квадратных матриц");
		if (precision < 0)
			throw new RuntimeException("Точность должна быть неотрицательной");
		
		// Начальный вектор - строчные суммы матрицы
		Matrix evector = new Matrix(rows(), 1);
		for (int i = 0; i < rows(); i++)
			for (int j = 0; j < columns(); j++)
				evector.set(evector.get(i, 0) + get(i, j), i, 0);
		evector.normalize();
		
		int iteration = 0;
		double abs = Double.MAX_VALUE;
		Matrix current = null;
		while (iteration++ < EIGENVECTOR_MAX_ITERATION && abs > precision) {
			current = multiplucation(this, evector, current).normalize();
			
			abs = 0;
			for (int i = 0; i < evector.rows(); i++)
				abs = Math.max(abs, Math.abs(current.get(i, 0) - evector.get(i, 0)));
			
			evector.copyFrom(current);
		}
		if (abs > precision)
			throw new MaxIterationException("Превышено допустимое число итераций: " + EIGENVECTOR_MAX_ITERATION);
		
		return evector;
	}
	
	public static Matrix addition(Matrix A, Matrix B) {
		return addition(A, B, null);
	}
	
	public static Matrix addition(Matrix A, Matrix B, Matrix Res) {
		return performElemOperation(A, ElemOperation.ADDITION, B, Res, true);
	}
	
	public static Matrix subtraction(Matrix A, Matrix B) {
		return addition(A, B, null);
	}
	
	public static Matrix subtraction(Matrix A, Matrix B, Matrix Res) {
		return performElemOperation(A, ElemOperation.SUBTRACTION, B, Res, true);
	}
	
	private static Matrix performElemOperation(Matrix A, ElemOperation operation, Matrix B, Matrix Res, boolean allocForRes) {
		if (A.rows() != B.rows() || A.columns() != B.columns() ||
			Res != null && (A.rows() != Res.rows() || A.columns() != Res.columns()))
			throw new IncompatibleDimensionsException();
		
		if (Res == null || allocForRes && (Res == A || Res == B))
			Res = init(new Matrix(), A.rows(), A.columns());
		
		for (int i = 0; i < Res.rows(); i++)
			for (int j = 0; j < Res.columns(); j++)
				operation.perform(A, B, Res, i, j);
		
		return Res;
	}
	
	private static interface ElemOperation {
		
		public static final ElemOperation ADDITION = new ElemOperation() {
			public void perform(Matrix A, Matrix B, Matrix Res, int row, int column) {
				Res.set(A.get(row, column) + B.get(row, column), row, column);
			}
		};
		
		public static final ElemOperation SUBTRACTION = new ElemOperation() {
			public void perform(Matrix A, Matrix B, Matrix Res, int row, int column) {
				Res.set(A.get(row, column) - B.get(row, column), row, column);
			}
		};
		
		public void perform(Matrix A, Matrix B, Matrix Res, int row, int column);
	}
	
	public static Matrix multiplucation(Matrix A, Matrix B) {
		return multiplucation(A, B, null);
	}
	
	public static Matrix multiplucation(Matrix A, Matrix B, Matrix Res) {
		if (A.columns() != B.rows() ||
			Res != null && (A.rows() != Res.rows() || B.columns() != Res.columns()))
			throw new IncompatibleDimensionsException();
		
		if (Res == null || Res == A || Res == B)
			Res = init(new Matrix(), A.rows(), B.columns());
		
		for (int i = 0; i < Res.rows(); i++)
			for (int j = 0; j < Res.columns(); j++)
				for (int k = 0; k < A.columns(); k++)
					Res.set(Res.get(i, j) + A.get(i, k) * B.get(k, j), i, j);
		
		return Res;
	}
}